Self-consistent calculation of the coupling constant in the Gross-Pitaevskii equation 
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A method is proposed for a self-consistent evaluation of the coupling constant in the Gross- 
Pitaevskii equation without involving a pseudopotential replacement. A renormalization of the cou- 
pling constant occurs due to medium effects and the trapping potential, e.g. in quasi- ID or quasi-2D 
systems. It is shown that a simplified version of the Hartree-Fock-Bogoliubov approximation leads 
to a variational problem for both the condensate and a two-body wave function describing the be- 
haviour of a pair of bosons in the Bose-Einstein condensate. The resulting coupled equations are 
free of unphysical divergences. Particular cases of this scheme that admit analytical estimations are 
considered and compared to the literature. In addition to the well-known cases of low-dimensional 
trapping, cross-over regimes can be studied. The values of the kinetic, interaction, external, and re- 
lease energies in low dimensions are also evaluated and contributions due to short-range correlations 
are found to be substantial. 

PACS numbers: 03.75.Hh, 05.30.Jp, 67.40.Db 



I. INTRODUCTION 

The Gross-Pitaevskii (GP) equation [1] is a power- 
ful tool for describing most of the physical properties of 
Bose-Einstein condensates of trapped alkali atoms (see 
reviews [2-5]). In the GP approach, the ground state en- 
ergy of a trapped dilute Bose gas of atoms of mass m is 
the functional 



E 



dr 



fi 2 
2m 



.F ext (r)H 2 + ||^ 



(1) 



of the order parameter <f> = <fi(r) = (W(r)), where ^(r) 
is the bosonic field operator and V eX t(r) is an external 
trapping potential. The coupling constant g in the GP 
functional (1) is intimately related to the density expan- 
sion of the energy of the homogeneous Bose gas. Indeed, 
in the homogeneous case Eq. (1) yields E/N = gn/2 
where n = N/V ~ \4>\ 2 is the 3D density. This expres- 
sion should be equal to the known first term in the density 
expansion E/N — 2i:h 2 an/m in three dimensions [6, 7]. 
With a being the 3D scattering length, the coupling con- 
stant g — Ai:h 2 a/m coincides with the zero momentum 
limit of the scattering amplitude, the two-body T-matrix, 
for two particles scattering in a vacuum. This standard 
approach based on the low-density expansion of the ho- 
mogeneous gas neglects the influence of inhomogeneous 
trapping potentials which may require a renormalization 
of the coupling constant. 

The situation is more complicated for two dimen- 
sional Bose gases, which can be regarded as the limit- 
ing case of a 3D gas with a highly inhomogeneous trap- 
ping potential. Kolomeisky et al. [8] proposed that the 
form (1) of the energy functional is still valid. How- 
ever, in this case the coupling constant becomes de- 
pendent on the local density. Indeed, Schick's result 
for the energy E/N = 2irH 2 n2D/[—m\n(n2DCi2 D )] of a 
dilute 2D Bose gas [9] leads to the coupling constant 
g = 47rfi 2 /(— ?7iln |0| 2 a2 D ) [8] and further corrections 



were derived in Refs. [10-13]. Here, n^o and a^o de- 
note the two dimensional density and scattering length, 
respectively. This generalization can be understood as 
a local density approximation, which yields consistent 
energy values in the homogeneous and inhomogeneous 
cases but does not reveal the nature of the additional 
non-linearity in the GP equation. Moreover, it is not 
clear what the nonlinearity should be in the crossover 
regimes from 2D to 3D and from ID to 3D. A mathe- 
matically rigorous justification of the GP functional [14] 
is of importance but hardly can help us in this situation. 

The purpose of the present paper is to show 
how a density-dependent renormalized coupling con- 
stant emerges naturally in a simplified Hartree-Fock- 
Bogoliubov (HFB) approximation starting from the bare 
interaction potential V(r). To this end, we derive gener- 
alized GP equations where the order parameter is coupled 
to the pair wave function of two bosons in the conden- 
sate. The latter has already been discussed in detail in 
Refs. [13, 15-17]. The generalized GP equations permit 
us to consider interaction potentials with a hard core di- 
rectly without the ^-function replacement by accurately 
accounting for the short-range spatial correlations of the 
particles. These correlations become essential in low di- 
mensions since the Born approximation for two scatter- 
ing particles fails at small momenta (see, e.g., Ref. [19], 
Sec. 45). We note that the correct treatment of the short- 
range correlations is also possible within the Jastrow[20] 
and the Faddeev-Yakubovsky [21] approaches. 

It is well known that the original HFB scheme leads to 
an artificial gap in the spectrum [22, 23]. Moreover, this 
scheme in conjunction with the (5-function replacement 
has an ultraviolet divergence [7, 24]. These problems 
then have to be cured by further approximations as classi- 
fied by Griffin [25] . Alternatively, complicated renormal- 
ization procedures [26, 27] or pseudopotentials [28, 29] 
have been suggested. In this paper we will discuss a 
novel approximation derived from the full HFB scheme 



where the use of the bare two-body potential provides an 
implicit renormalization, and the ultraviolet divergences 
are avoided. We will discuss the excitation spectrum and 
show that it is gapless in a reasonable approximation. 

Low dimensional Bosc systems are not only of gen- 
eral theoretical interest but also find the attention of 
current experimental exploration [30-33]. In these ex- 
periments, the low dimensional condensate is realized 
by highly anisotropic 3D trapping potentials when the 
single-particle energy-level spacing in the tightly confined 
dimensions exceeds the interaction energy between atoms 
hw PtZ > [i. Here the frequencies are associated with the 
axially symmetric harmonic potential, and /i is the part of 
the chemical potential coming from the particle interac- 
tion, which is of order of the mean interaction energy per 
particle. This criterion takes the form l p ^ z < £ in terms 
of the coherence length £ = H/ \fjixn [34] and the radial 
(axial) harmonic oscillator length l p z = y / /j/(mw p . z ). 

Theoretically, the 2D regime l z <C £, l p 3> £ was in- 
vestigated in detail by Petrov, Holzmann, and Shlyap- 
nikov [35] . The coupling constant was assumed to be the 
T-matrix of two particles scattering in the harmonic trap 
with l p — oo at the energy of relative motion E = 2/i. 
An additional nonlinearity is introduced, since the local 
value of \x depends on the density and the coupling con- 
stant in a self-consistent way. In this regime, the motion 
of particles is confined in z-dircction to zero-point os- 
cillation. This implies that the order parameter can be 
represented in the form <f>(x,y,z) — <j>o(z)<p(x,y), where 
<fio(z) is the ground state of the ID harmonic oscillator 
and 4>{x, y) is governed by the two dimensional GP equa- 
tion resulting from the functional (1) in two dimensions. 
So, in this regime, the behaviour of the condensate in x-y 
plane is the same as in the "pure" 2D case with the 2D 
scattering length written in terms of the length l z of the 
tight confinement [35]. 

An improved many body T-matrix theory was devel- 
oped by Stoof and coworkers [36] in order to describe not 
only the homogeneous low-dimensional Bose gases but 
also the crossover from 3D to lower dimensions. The cou- 
pling constant in the inhomogeneous case is represented 
by the local value of the T-matrix, which depends on the 
local value of the chemical potential. The local T-matrix 
approximation was also used in Rcf. [37] and the micro- 
scopic approach of Ref. [38]. Thus, one can say, slightly 
simplifying the situation, that the common method of 
evaluating of the coupling constant in the above works is 
to determine first the T-matrix from the corresponding 
two-body Schrodingcr equation supposing that the mo- 
tion of the particles is infinite in some directions, and 
after that to replace the coupling constant by the local 
value of the T-matrix. In this paper we offer a method 
beyond the local T-matrix approximation. The coupling 
constant is determined self-consistently for a given 3D 
geometry from a unified variational scheme. As a re- 
sult, we obtain a non-local term in the energy functional, 
which can be of practical importance if the external po- 
tential varies on the scale of the interaction potential. 



This may be realized, e.g., for condensates of loosely- 
bound molecules in tight or strongly oscillating potentials 
like optical lattices. We expect experiments to enter this 
regime in the near future as both atomic condensates in 
optical lattices [31, 39] and molecular condensates [40-42] 
are currently under intense experimental investigation. 

As a starting point of our approach we assume that 
we have a Bose-Einstein condensate or quasi-condensate 
with a well-defined order parameter. Long-range fluctu- 
ations of the phase, which become important for many 
physical properties in low-dimensional traps [43] , are be- 
yond the scope of our scheme. They can be studied by 
means of the approaches of Refs. [35, 36, 44, 45]. Also 
the strongly- interacting fermionized regime of the Tonks- 
Girardeau gas [32], which was studied in Ref. [46], can- 
not be described with the methods of this paper. How- 
ever, we note that our scheme, within its accuracy, is 
simple and physically transparent and able to reproduce 
not only the value of the coupling constant in ID [47] 
and 2D [35] regimes but also to describe the 3D-2D and 
3D-1D crossovers. Furthermore, it allows us to calculate 
directly the correct values of the kinetic and interaction 
energies of bosons in the trap, which are not given by the 
first and the third terms, respectively, in the GP func- 
tional (1) [48]. 

The paper is organized as follows. In Sec. II we derive 
the generalized GP functional and corresponding equa- 
tions from a simplified HFB approximation. In Sec. Ill, 
a few specific cases are considered that admit analytical 
estimations, including the homogeneous and inhomoge- 
neous Bose gases in low dimensions. In Sec. IV we cal- 
culate the values of various contributions in the energy 
In particular, the release energy of the low dimensional 
gases is estimated. In Sec. V we derive a useful virial 
theorem and a relation between the chemical potential 
and different parts of the energy functional. The eigen- 
functions of the two-body density matrix and a relation 
between the normal and anomalous averages are obtained 
within the HFB approximation in Appendices A and B, 
respectively. 



II. GENERALIZED GROSS-PITAEVSKII 
EQUATIONS 

A. Failure of standard GP approach 

In the standard approach [2-5], the equilibrium value 
of the order parameter (f> is determined by minimization 
of the GP functional (1) with the constraint of particle- 
number conservation S(E — /i' N) / 6<f>* (r) = 0. Here N ~ 
Nq is the number of particles and the chemical potential 
[t! appears as a Lagrange multiplier. Introducing for later 
convenience ji = fj,' — Eq as the chemical potential due 
to interaction where Eq is the ground state energy of a 
non-interacting particle, we arrive at 



(E + /z)</> : 



2m 



V 2 + y cxt (r)0 + ff |0| 2 0. (2) 



The simplicity of this derivation is based on the simple 
form of the GP energy functional (1), where the effects 
of the binary inter-particle interactions has been reduced 
to a single parameter given by the coupling constant. 
In order to determine the constant self-consistently, it 
should be examined carefully how the interaction term 
(g/2)|0| 4 appears in the GP functional (1). 

In a general many-body system with binary inter- 
actions, the expectation value of the interaction en- 
ergy is a functional of the two-body density matrix 
(W (xi)W (x 2 )^> {x' 2 )^ (x'i)) ■ For a pairwise interaction 
potential V(a;i,a:2) = V"(ri — r2,<7i,<72) we thus obtain 
[19, 49] 



Eint = (-'^2v(x i ,x : j)\ = - dx 1 dx 2 V(x 1 ,x 2 ) 



x(^(x 1 )^ t (x 2 )^(x 2 )*(x 1 )), 



(3) 



where x = (r, a) stands for the coordinate and spin or 
sort indices of a particle, respectively and J dx ■ ■ ■ = 
5^ CT J dr • • • . The kinetic energy and the energy of in- 
teraction with an external field are determined by the 
one-body matrix (^'(x)^!(x')) 



Ev\ 



kin 



-Efixt, — 



E 



2m 



— / dx Vl(&tf)ir(x)) 

Ira 



(4) 



J2 V ***( X *))= [dxV e ^(x)(&(x)i!(x)).(5) 



The behaviour of the one- and two-body matrices is 
easy to understand in the dilute limit, when the con- 
densate depletion (N — No)/N is small. We note that 
the number of bosons in the condensate iVo is defined 
as the macroscopic eigenvalue of the one-body density 
matrix (&(x)${x')), that is J dx' {&(x')$(x))<j>o(x') = 
No<po(x). The field operator can be expanded in the com- 
plete orthonormal set of eigenfunctions of the one-body 
matrix 'I'(x) = ao<po(x) + "^2 „o„^„(a;), where the sum 
Y^, v means X^o an< ^ /dx |</> y (x)| 2 = 1. Appearance 
of the Bose-Einstein condensate implies the macroscopic 
occupation of No, i.e. the ratio Nq/N remains finite in 
the thermodynamic limit. Following Bogoliubov [6, 50] 
we now replace the condensate operators by c-numbers 
Oq = ao — VNq and represent the Bose field operator in 
the form *(x) = <j>(x)+&{x) [51]. Here <j>(x) = y/No<l>o(x) 
is the c- number part, and d(x) = ^ v a v <j) v {x) the 
operator part, for which we have ('I'(x)) = <p(x) and 
(i?(x)} = 0. Thus, the order parameter is nothing else 
but the non-normalized eigenfunction of the one-body 
density matrix. 

In the original approach of Gross and Pitaevskii, the 
simplest mean-field approximation is used when the op- 
erator part is completely neglected: Vt'(x) ~ <fr(x) and 
$!>(x) — 4>*{x). Assuming additionally that the order 



parameter does not change significantly at the distances 
of order of the radius R c of the interaction potential, we 
obtain the GP energy functional (1) for spinless bosons 
from Eqs. (3-5) with the coupling constant 



g B = / drV(r). 



(6) 



This coupling constant can be identified with the two- 
body scattering amplitude at zero momentum in the 
Born approximation. The validity of the GP approach 
with the coupling constant g B of Eq. (6) is certainly 
linked to the validity condition of the Born approxima- 
tion at zero momentum that the potential V(r) be small 
and integrable. 

Of the two assumptions mentioned above, namely the 
slow spatial variation of the order parameter and the va- 
lidity of the Born approximation, the former is usually 
fulfilled as the healing length £ = H/^/JIm as a lower 
bound of the length scale of the order parameter [3, 4] 
is usually much larger than the effective range of the in- 
teraction. The latter assumption, however, is clearly not 
fulfilled for experiments with cold atomic gases as their 
interactions are of the hard-core type. 

The validity of the GP approach can be extended 
to such systems by an argument attributed to Lan- 
dau [6, 16]. He noted that at extremely low energies, 
as predominant in the dilute-gas BEC, the scattering 
properties are completely determined by only one single 
parameter, which is the 3D s-wave scattering length a. 
This allows us to replace the Born approximation for the 
scattering amplitude g B by its exact value g = 4irH 2 a/m, 
which can be found from the two-body Schrodinger equa- 
tion even for hard-core potentials. This indirect argu- 
ment, however, cannot be used in one or two dimensions 
as there is no such simple relation between the integral 
in Eq. (6) and the scattering amplitude as we have in 
three dimensions. Furthermore the Born series for the 
scattering amplitude diverges for small momenta in two 
dimensions and below (see, e.g., Ref. [19], Sec. 45). 



B. Pair wavefunction in a medium 

The above described deficiencies of the naive GP ap- 
proach may be remedied by accounting for the two- 
particle scattering processes explicitely. Within the HFB 
scheme this is possible through certain correlations intro- 
duced by the fluctuation operators i). In order to see the 
relation between the two-particle scattering process and 
the correlation functions mentioned above it is useful to 
introduce the concept of a two-body or pair wave func- 
tion in the medium of other particles [15, 50]. The pair 
wave functions in the medium of the many-body system 
arc defined as eigenfunctions of the two-body density ma- 
trix, as discussed in detail in Appendix A. They should 
be distinguished from the two-body wave functions in the 
vacuum, which relate to a system of two particles. For 
the latter we will use the superscript (°) . 



Let us suppose that we know the exact eigenfunctions 
of the two-body density matrix. Then we can expect for 
the dilute gas, where low-momentum two-body processes 
dominate the behaviour of the system, that the pair wave 
functions in the medium should be very close [15] to the 
two-body wave functions in the vacuum, which are the 
solutions of the two-body Schrodingcr equation. This 
physical assumption was used to obtain the density ex- 
pansions for the 3D [16, 17] and 2D [13] homogeneous 
Bose gases in a very simple manner. However, various 
approximations in the many-body theory can break this 
relation. 



C. A simplified HFB scheme 

Within the HFB approximation for the homogeneous 
Bose gas, all eigenfunctions of the two-body density ma- 
trix except for one are plane waves and are thus treated in 
the Born approximation as shown in Appendix A. This 
is an obvious drawback of the HFB scheme. It turns out 
that the two-body wave function that is not a plane wave 
is proportional to the anomalous average 



if(x u x 2 )= ($(xi)$(x 2 )) 



(7) 



and corresponds to the macroscopic eigenvalue Nq(Nq — 
1) in the limit of large N. It is the pair wave function 
that describes the two-particle scattering process in the 
medium of the Bose-Einstein condensate [15, 52]. Thus 
we can go beyond the Born approximation in the frame- 
work of the HFB method by keeping only the contri- 
bution of the anomalous average (\1/ (xi)$ (#2)) and ne- 
glecting the contribution of the other wave functions in 
the two-body density matrix. Due to small condensate 
depletion (N — N )/N -C 1, one can expect that the con- 
tribution of only this eigenfunction will be sufficient for 
obtaining the coupling constant in the GP equation. In 
this simplified version of the HFB approximation we set 

(^(X 1 )¥(X2)^(X' 2 )^>(X' 1 )) ~ V *(x 1 ,X2Mx' 1 ,x' 2 ). (8) 

Extracting the c-number part of the field operator, the 
anomalous averages can be rewritten as 



ip(xi,x 2 ) = (j){xi)<f)(x 2 ) +ip(xx,x 2 ), 



(9) 



where we introduced the notation ip(xi,x 2 ) = 
( , &(xi)'d(x 2 )) for the anomalous two-boson correlation 
function associated with the scattering part of the 
two-body wave function. The functions (p(xi,x 2 ) and 
ip(xi,x 2 ) are symmetric with respect to permutation of 
X\ and x 2 due to the commutation relations for the Bose 
field operators. 

For the one-body density matrix we find 
(^(x)i'(x')) = <j)*{x)4>{x') + (tft(x)tf(x')). We 
note that the normal ($' (x)'0(x')} and anomalous 
('d(x)'d(x')} averages are not independent quantities as 
discussed in Appendix B and Refs. [16, 17]. Within the 



Hartree-Fock-Bogoliubov scheme, the relations between 
them appear as a specific property of the HFB ground 
state (the quasiparticle vacuum) and do not contain 
parameters of the Hamiltonian in explicit form. We will 
use the approximate relation (B20), which leads to 

(&(x)$(a/)) =0*( x ) ( j ) (x')+ [ dx 2 i>*(x,x 2 )i;(x 2 ,x'). 



(10) 
With the help of Eqs. (8), (9), and (10) we rewrite 
Eqs. (3), (4), and (5) in terms of the anomalous aver- 
ages 

E int = - dx 1 dx 2 V(x 1 ,x 2 )\(p(x 1 ,x 2 )\ 2 , (11) 



E cxt = ^ / dx 1 dx 2 [V cxt (x 1 ) + V cxt (x 2 )]\tp(x 1 ,x 2 )\' 



+ / dariKxtfriMari)!' 



1 



(12) 



^kin = ■= J dxidx 2 ip*(x 1 ,x 2 )(T 1 +T 2 )ip(xi,x 2 ) 



+ / dxi 4>*{xi)Ti(f>(xi), 



(13) 



where Tj = — h 2 V 2 /(2m) and j = 1, 2. The total number 
of particles is related directly to the one-body matrix: 
N = fdx{&(x)$(x)). WithEq. (10) we find' 

N = f d»i |0(xi)| 2 + f dx 1 dx 2 |V(zi, x 2 )\ 2 . (14) 

The current approximations are useful for a variational 
scheme where the functions <j)(xi) and ip(xi,x 2 ) are de- 
termined by minimization of the total energy with the 
constraint N = const. Using the Lagrange method, 
we obtain the conditions 8E/84>(xi) = 8E/8<j>*(x\) — 
SE/8ip(xi,x 2 ) — 8E/8ip*(x\,x 2 ) — for the energy func- 
tional 

E[{<f>, V>}, n'] = E kin + E cxt + E int - n'{N - AT), (15) 

given by Eqs. (11)-(14). Here // = /j, + Eq is the chemical 
potential and Af = (N) is the average number of parti- 
cles, i.e. the l.h.s. of Eq. (14) at the equilibrium values 
of the functions <p and -0 corresponding to the minimum 
(ground state) of the functional (15). Note that the varia- 
tion 8ip(xi, x 2 ) is symmetric under the permutation of x\ 
and £2, so, the equation J dxidx 2 g(xi,x 2 )8ip(xi,x 2 ) = 
leads to g(xi,x 2 ) + g{x 2 ,xi) = for arbitrary functions 
g(xi,x 2 ). 

This variational procedure yields the following system 
of equations for the one- and two-body functions (j>(xi) 
and ip(xi,x 2 ), respectively, 



Ci<j)(xi) + / dx 2 (/)*(x 2 )V(xi,x 2 )ip(xi,x 2 ) = 0,(16) 
(Ci+C 2 )i/)(xi,x 2 ) + V(xi,x 2 )ip(xi,x 2 ) = 0,(17) 
where the operators C\ and C 2 stand for 

Cj = -h 2 V 2 /{2m) -n-E + V ext ( Xj ), j = 1, 2, 



and <j), ip and ip are simply related by Eq. (9). Due to 
this relation, Eq. (16) is nonlinear with respect to <j> and 
can be associated with the GP equation. Equation (17) 
is the analogue of the two-particle Schrodinger equation 
and is linear with respect to ip, though not uniform. The 
obtained system of two equations allows us to determine 
the coupling constant self-consistently. A specific fea- 
ture of Eq. (16) is the non-local nature of the last term, 
which can play a role when the radius of the interact- 
ing potential becomes of the order of the characteristic 
length of the anisotropic trapping potential in some di- 
rection, say, R c ~ l z <C £, or if the trapping potential 
has a strongly oscillating contribution with the scale of 
the order of R c . At the same time, Eqs. (16) and (17) 
indeed reduce to the GP equation with the 3D coupling 
constant g — 4irh 2 a/m in the limit R e -C £ -C l x , l y , l z 
as will be shown in Sec. IIIB 1. 

When the external potential becomes independent of 
some coordinate, say z, particles can move freely in z- 
direction and we should impose the boundary conditions 
that follow from Bogoliubov's principle of correlation 
weakening [50]: 

$(x)${x')) ~ ($(x)){${x')) when \z - z'\ -> 00. 

Physically, this implies that the function ip(xi,x 2 ) = 
( , d(xi)'d(x 2 )) vanishes at the spatial distances of order 
of the coherence (healing) length, |ri — r^l > £. 

A time-dependent generalization of Eqs. (16) and (17) 
can in principle be derived from the equations of mo- 
tion of the field operators. Here, however, we will not 
elaborate the full derivation but instead present a simple 
argument that leads to a useful time-dependent scheme. 
In the case of a time-independent Hamiltonian, it can be 
easily seen that the GP order parameter depends on time 
as 

4>(x,t) = (N-l\$(x,t)\N) 

= 4>(x)cxp[-i(E N - E N -i)t/H] 
— (f>(x)exp[—ifi't/h]. 

Here, \N) and En arc the ground state and en- 
ergy of N bosons, respectively. By analogy, we find 
ip(x\,x 2 ,t) — ip(xi, X2) exp[— i2fi't/H] and (p(xi,x 2 ,t) = 
(p(xi, X2) exp[—i2/j,'t/fi}. We now argue that the chem- 
ical potential in Eqs. (16) and (17) arises due to time 
derivatives, which leads to the obvious generalization 



d 

ih—ip(xi,X2,t) 



where we denote 



Hi(j)(xi,t) + E n i(xi,t), 



(18) 



[Hi + H 2 + V{x\,x 2 )]ip(x\,x 2 ,t) 

+(f)(xi,t)E nl (x 2 ,t) 
+<Kx2,t)E n i(xi,t), (19) 



H, 



2m 



+ V vxt (x j ,t), J = 1,2, 



E n \(x,t) = / dy 4>*(y,t)V(x,y)ip(x,y,t). 



The functions = (&(x,t)} and <p = (\I>(xi,£)\I>(x2,t)) 
are normalized as J dx \</>(x, t)\ 2 = Nq and 
J dxidx 2 \ip(xi,x 2} t)\ 2 — N (N — 1) ~ Nq, respectively. 
The time-dependent generalized GP equations (18-19) 
become the ordinary one- and two-body Schrodinger 
equations, respectively, in the limit £ 3> l x ,ly,lz when 
we can neglect all the nonlinear terms, which are 
responsible for many-body effects. Therefore they are 
of slightly more general validity than the stationary 
equations (16-17), which imply a large particle number, 
since E N — E N - 2 — 2(En — En-i) — 2// is valid 
only for N ^> 1. We notice that the time-dependent 
equations similar to that of (18) and (19) were derived 
in papwer [67] by the method of noncommutativc 
cumulants. 



D. Properties and limits of validity 

The time-dependent Equations (18-19) give access to 
the elementary excitation spectrum. At the moment 
we cannot prove the gaplessness of the spectrum in the 
most general case, but we can solve for the excitation 
energies approximately. With the ansatz ^j(ri,r2,i) = 
4>(ri,t)<fi(r 2 ,t)[l + i/j(r)/no], we obtain the Bogoliubov 
excitation energy lou = y/T% + 2nolI(k)Tk with the k- 
dependent scattering amplitude U(k) (for the notations 
see Sec. Ill A). This form of the spectrum for a singular 
two-particle interaction was proposed without derivation 
by Bogoliubov in Ref. [6]. For small k we can replace 
U(k = 0) = A-Kh 2 a/m and obtain the usual (gapless) 
Bogoliubov dispersion. The additional features in the 
obtained spectrum at medium and high energies reflect 
the structure of the interaction potential neglected in the 
standard GP approach and present a clear advantage of 
our extended scheme. As a consequence, we can expect 
that Levinson's theorem for quasi-particle scattering [53] 
will be modified. 

Let us discuss limits of validity of the generalized GP 
equations (16) and (17). First, we imply that the Bose- 
Einstcin condensate (or quasi-condensate in low dimen- 
sions, see Sec. Ill B) is developed strongly. This means 
that ro ^C £, where r is an average distance between 
bosons [3, 4]. Second, the above derivation can be ap- 
plied only to the short-range interaction potentials that 
decrease at least as fast as V(r) ~ \/r e+D for r — ► 00, 
where D is dimension and e > [16, 17]. For a long- 
range interaction like Coulomb repulsion, the approx- 
imation (8) works badly. Third, the approximations 
(8) and (B20) arc insufficient to describe the long-range 
behaviour of the normal (i?t (xi) , d(x 2 )) and anomalous 
(t?(xi)i?(x2)) correlation functions, which are governed 
by Bogoliubov's "1/g 2 " theorem [50, 54, 55]. According 
to this theorem, the above correlation functions should 
decay as l/|ri — r2 1 2 when |ri — r2 1 > £ at zero tempera- 
ture if the Bose-Einstein condensate exists. Our scheme 
gives l/|i"i — r2 1 decay, as we show in Sec. Ill A. How- 
ever, we stress that the long-range behaviour is not needed 



for obtaining the coupling constant, since the integral 
in Eq. (16) contains the anomalous correlation function 
multiplied by the short-range potential V{x\,X2) with 
the characteristic radius R c < £. Since the developed 
scheme describes well only the short-range behaviour of 
tp(x\ , x-i) for |ri — V2 | < £, the integration in the last term 
of Eq. (14) should be restricted to this region 



N= / dx 1 \0(x 1 )\ 2 + 



dxidiE2 \ip{x 1 ,x 2 )\ 2 , 

(20) 
otherwise we obtain formally divergent term. This modi- 
fication of the original scheme, however, does not change 
the working equations (16 - 19) in the region |ri — r2| ^ £, 
which is of sole interest for our purposes. We stress that 
Eq. (20) is really needed only when minimizing the en- 
ergy functional (15) directly. Furthermore, if we obtain 
the solutions of Eqs. (16) and (17) as functions of the 
chemical potential in the grand canonical ensemble, then 
the condition (14) or (20) can be employed without the 
second term at all in order to rewrite the answer in terms 
of the total number of particles in the canonical ensemble, 
owing to small condensate depletion. 

Note that the standard HFB approximation can be 
obtained by using the variational scheme if, first, one 
substitutes Eqs. (3)-(5) into the energy functional (15), 
second, employ the restrictions Eqs. (B18) and (B19), 
and third, retain all additional terms missing in Eq. (8), 
where the three- and four-boson averages of d and if' 
ought to be evaluated by means of the Wick's theorem 
and, consequently, the three-boson averages vanish. 



III. EXAMPLES 

In this section we restrict ourselves to spinless bosons 
with an isotropic short-range interaction V = V(r), 
where r — |ri — i"2 1 - Even after this simplification, the 
solution of the generalized GP equations (16) and (17) 
remains a rather complex problem. Nevertheless, in a 
number of specific limiting cases we are able to obtain 
analytic results. 



A. The homogeneous case 

Let us investigate Eqs. (16) and (17) in three and two 
dimensions for the homogeneous Bosc gas. In the homo- 
geneous case V cxt — 0, hence we have ip — ip(r), E = 0, 
and Eq. (16) gives the trivial solution cj> = ^JfMy = const. 
In this subsection, we use the common notation no for 
the condensate density in both 2D and 3D cases. Thus, 
Eqs. (16) and (17) read 

" = / drV(r)[n + ip(r)], 

h 2 
2ixip{r) = V 2 ^{r) + V{r)[n + ip{r)}, 



and ip(r) — ► for r — > oo in accordance with Bogoliubov's 
principle of correlation weakening. Taking the Fourier 
transformation of the last equation, we obtain 



n 



= n Q U(0), (21) 

= -P.P. ^ r, (22) 

where we denote U(k) = j dr V(r)e _ r [l + t/j(r)/no], 
Tk = H 2 k 2 /(2m), and the symbol P.P. stands for the prin- 
ciple value of the associated integral. The latter appears 
as a natural regularization for the singular denominator 
in the r.h.s of Eq. (22) and implies that the scattering 
part of the two-body wave function ip(k) is real and cor- 
responds to a standing wave. Note that another regular- 
ization, such as the standard replacement k — > k ± is, 
leads to the same results in the leading order at small 
densities. Within the more accurate method [16, 17], we 
obtain the same equation as (22) but with the Bogoli- 
ubov denominator 2^/T^ + 2noU(k)Tk. The latter pro- 
vides the correct values of both the short- and long-range 
behaviour of the correlator ip(r) — (#(r)#(0)) [which is 
the Fourier transform of tp(k)], while Eq. (22) provides 
only the short-range behaviour. Indeed, in the 3D case 
we have tp(r) <~ cos(\/2V/£)/r at r > £ (see below) but 
not ip{r) ~ 1/r 2 as it should be. 

Equation (22) can be rewritten in the Lippmann- 
Schwinger form with the help of the Fourier transforma- 
tion. By using the familiar property of Fourier transfor- 
mation /dke 4kr 5 (k)/(k)/(2 7 r)' D = / dr'/(r')5(r - r') 
(here D is the dimension), we obtain the equation for 
ip(r) = n + ip(r) 



<p(r) - n + / dr' V(r')<p(r')G(\r - r'|), 



where the Green function is introduced 



G(r) = -P.P. 



dk 



oik-r 



(2Tr) D 2{T k -^)- 



(23) 



(24) 



In the dilute limit, when the average distance between 
particles is much less than the coherence length, the wave 
function ip(r)/no, describing the behaviour of two parti- 
cles in the condensate, should be proportional [15, 56] to 
the s-wavc function ip(°>(r), which corresponds to relative 
motion of two particles with zero momentum and obeys 
the two-body Schrodingcr equation in the center-of-mass 
system 



-(fi7m)VV V) +V(r)<p w (r) = 



(°)M = 



(25) 



In the 3D case, the coefficient of proportionality is 
equal to unity [6] in the leading order with respect to 
the density, provided the following boundary conditions 
arc imposed: first \if^ (r)\ < oo at r = and second, 

</>3rj(r) — > 1 — a/r for r — > oo. In the developed for- 
malism, this can be easily inferred from the obtained 
equation (23). Indeed, direct integration in Eq. (24) 



gives Gsuir) = — TOCOs(v / 2? , /£)/(47r?i 2 r), and, hence, 
G3d(j) — — m/(4irH 2 r) when r < £. Thus we have 
<p(r) ~ no(p 3 jy(r) within this region, and integration of 
Eq. (25) yields U(0) = 47rfi 2 a/m. For the dilute gas we 
have also Uq ~ n, and Eq. (21) leads to the familiar ex- 
pression for the chemical potential /j, ~ 4irH 2 na/m. 

In the 2D case, the low-energy behaviour of the 2D 
Green's function (24) is easily calculated: G 2 d(t") — 
m/(27r?i 2 ) ln[e 7 r/(V2^)] when r < £. Then it is not diffi- 
cult to see from Eq. (23) that, first, <p(r)/no obeys the 2D 
Schrodinger equation (25), and, second, its asymptotics 
for r — > oo is 

ip(r)/no -+ 1 + ln[e 7 r/(V20]™C/(0)/(2vrfi 2 ). (26) 

Hence, due to linearity of Eq. (25), the solution for <p(r) 
should be proportional to the wavefunction ^dM that 
obeys the 2D Schrodinger equation (25) with the follow- 
ing boundary conditions: first |<p 2r j(r)| < oo at r = sec- 
ond, <^2D( r ) ~~* m ( r / a 2D) for r — ► oo. The latter equation 
can be considered as the definition of the 2D scattering 
length [57]. Note that in the case of hard disks, <z 2 d co- 
incides with the radius of the disks. It is convenient to 
introduce the dimensionless parameter u by the relation 
U(0) = 4:TTh 2 u/m, such that u is the dimensionless scat- 
tering amplitude for two bosons in a medium of other 
bosons. By comparing the asymptotics (26) with that of 
(/5 2E J(r), we derive 



tp(r) 
ln(a 2 D/0 



l/(2w)+ln(e 7 /\/2). 



(27) 



With the help of Eq. (21) and the definition of £ (see 
above), the relation (27) becomes a self-consistent equa- 
tion for u 



1/u + \nu = — \n(n2Da 2 i)2TT) — 2j, 



(28) 



where we neglect the condensate depletion in the lead- 
ing order, putting no ~ n 2 D- By means of the latter 
approximation, the expression (21) takes the form 



/i = 4ttH n2r>u/m. 



(29) 



Thus, the 2D chemical potential is given by Eqs. (28) and 
(29), which lead to the density expansion 



47rfi 2 n 2 D 



1 



1 



it) 



ln(n 2D a 2 



2DJ 



In 2 (n 2D a 2D ) 



x In 



ln(n 2 D»2E)) 



(30) 



Equations (28) and (29) are in agreement with the results 
of Refs. [10-12] and with the more accurate scheme of 
Ref. [13], which yields the correction for the chemical 
potential 



/i = (AnH n2T>/m)(u- 



(31) 



Here, u is given by the more exact relation 

l/u + lnw = — ln(n 2 Da2 D 7r) — 27, (32) 

where 7 = 0.5772 ... is the Euler constant. By means of 
this relation, one can rewrite Eq. (31) in terms of the gas 
parameter n 2 DQ 2 rj and obtain three more terms in the 
expansion (30). Note that Eq. (32) differs from Eq. (28) 
by a numerical factor under the logarithm, which is es- 
sential only for obtaining these additional terms but not 
the terms given by relation (30). 

B. The inhomogeneous case 

1. The Gross-Pitaevskii regime 

First of all, we should verify that the equations ob- 
tained in Sec. II lead to the standard GP scheme in the 
case R c -C £ <C I, where / is the characteristic length of 
an isotropic trap. In this regime, one can expect that 
the pair wave function </?(ri, r 2 ) is very close to that ob- 
tained in the homogeneous case, with the difference that 
the density is spatially dependent now. So, we put by 
definition ip(r 1 ,r 2 ) = </>(ri)</>(r 2 )^(ri,r 2 ) and^(ri,r 2 ) = 
0(ri)0(r 2 )^(ri,r 2 ), and, hence, p(n,r 2 ) = l+tp(r 1 ,r 2 ) 
by Eq. (9). Substituting those expressions into Eqs. (16) 
and (17) yields 



H 2 

— V 2 -M-^o + ^ext(ri) 
2m 



ri 



■0(n)|0(n)| s 



2m 



x /'dr 2 V(|p 1 -r 2 |)^(p 1 ,p 2 ) = ) (33) 

(V 2 + V|)^(ri,r 2 ) + V(ri - r 2 )£(n,r 2 ) 
= [/(ri) + /(r 2 )]^(n,r 2 ), (34) 



where we use the condition R e <C £ in the first equation 
and introduce the notation 



/W 



dr'\(f>(r')\ 2 V(r-r')lp(r,r') 



ft 2 V P 0(r) 



•v r , 

(35) 



with the last term being a differential operator. Since 
<^(ri,r 2 ) ~ ^(pi)<^(p 2 ) at the distances of order of the 
correlation length, we have £>(ri,r 2 ) ~ 1 at these dis- 
tances. Consequently, the l.h.s. of Eq. (34) remains finite 
when the density tends to zero, while the r.h.s. becomes 
small. Indeed, the first term of Eq. (35) is of order of 
h 2 an/m. The second term is less than h 2 /(m£, 2 ) because 
the characteristic scale of the order parameter cannot be 
smaller than £ in the case £ <C / and the same applies 
to tp. Hence, in the leading order we can completely ne- 
glect the r.h.s. of Eq. (34), which leads to the standard 
Schrodinger equation (25) for yp. Thus, we come to the 
approximation 



p(ri,r 2 ) ~ </>(ri)</>(r 2 )^(r 



(36) 



Using the well-known relation for the 3D scattering 
length 



47rfi 2 a/m= / d 3 r V(r)ip^(r) 



(0), 



(37) 



we can rewrite Eq. (33) in the standard GP form with 
the coupling constant g — 4irh 2 a/m. Note that, never- 
theless, the equilibrium value of the energy (15) differs 
from that of the GP value (1) by the terms arising from 
the condensate depletion because the second term in the 
r.h.s. of Eq. (10) is not equal to zero. We will discuss 
these corrections to the energy in Sees. IV and V. 



2. 2D regime 

Here we consider the Bose gas confined only in z- 
direction by the trapping potential V cx t — V cx t(z). The 
system is homogeneous in the x-y plane and assumed to 
be infinitely large. Physically this means that the x-y 
size of the system is much larger than the characteristic 
radius of the trapping potential l z = yjh/(mtjj z ). The 
order parameter (f> now becomes independent of x and y, 
and the two-body function depends on the relative dis- 
tance p — \p\ — P2I between points p\ — {x\,y\) and 
P2 = (£2,2/2), so ip(r 1 ,r 2 ) = (f(zi,z 2 ,p). The 2D regime 
is provided by the condition l z <C £. Moreover, the con- 
dition R c <C £ is fulfilled in most experiments. As was 
discussed in Sec. I the density profile is then governed 
by the ground state solution <po( z ) of the one-particle 
Schrodingcr equation 

[ - h 2 V 2 /(2m) -E a + V ext {z)]Mz) = 0, 

because the second term in Eq. (16) can be treated as 
a small correction. Thus, we can put in the leading or- 
der <j)(z) ~ y / n2D</ ) o( z ); 4>o{ z ) is normalized to unity. By 
analogy with standard perturbation theory, the chemi- 
cal potential, as the first correction to Eq, can be found 
with the unperturbed eigenfunction cf>Q. So, multiplying 
Eq. (16) by (j>o(zi) and integrating by z\ yield 



M = ™2D / dpU{p), 
where by definition 



(38) 



U(p) = / dz l dz 2 V(^p 2 + {z l -z 2 ) 2 ) 

Xip(zi,z 2 , p)(f>o(zi)4>o(z 2 )/n 2 B- (39) 

In the same manner, one can multiply Eq. (17) by 
4>o(zi)(j)o(z2) and carry out the integration by z\ and z 2 , 
which results in the equation 



2[n 2 V 2 p /{2m)+p]i,(p) = U{p) 
for the function 



(40) 



V>(/>) = dzidz 2 -il>(z ll z 2 ,p)4> (z 1 )4> (z 2 )/n 2 r>- (41) 



Thus, we arrive at the system of equations (38) and (40), 
which coincides with that of (21) and (22) in homoge- 
neous case if we put U(k) = J dp U(p)e~ t]i ' p and perform 
the Fourier transformation of Eq. (40). By the same 
method as in Sec. Ill A, we obtain the asymptotics for 
sufficiently large p [physically, for R c <C p <C £, when 
only the first term dominates in Eq. (40)] 



<p{p) ~ 1 + ln[e-'p/(V20}mp/{2irH 2 n 2D ), 
where by definition 



<P(P) 



dz 1 dz 2 ip(z 1 , z 2 , p)(j) (z 1 )(j) (z 2 )/n 2I) 



(42) 

l+^(p). 

(43) 
The latter relation is due to Eqs. (9) and (41). 

In order to obtain the chemical potential in terms of the 
3D scattering length a and the length l z of the trapping 
potential, we use the following approximation [58] 



,(0) 



ip(zi,z 2 ,p) = CV3 E j(r)n 2 D<M2 ; i)<M* z 2) 



(44) 



in the region r <C l z <C £, where r = \J p 2 + (z\ — z 2 ) 2 , 
and tp 3 -^ (r) denotes the 3D solution of the Schrodingcr 
equation (25) with asymptotics for r 3> R c 



^iSW 



1-o/r 



(45) 



Here the crucial point is that the constant C =/= 1, which 
determines the 2D behaviour of the system. If we substi- 
tute Eq. (44) into Eq. (43) and take the integral, we arrive 
at a new expression for <p(p). This should be expanded 
with respect to the dimensionlcss variable p/l z and com- 
pared with Eq. (42). Since the main contribution in that 
integral comes from the asymptotics (45), one can use 

it instead of the function tf^T) ( r ) itself. By performing 
this procedure for the harmonic trapping potential with 
4>q(z) = exp[-z 2 /(2/ 2 )]/ v /Z^0F, we have 



(p(p)^C + 



2Ca 

lzV2^ 



ln[e^ 2 p/(2V2l z )}. 



Comparing this relation with Eq. (42) yields 
C = y2wl z u/a, 



(46) 



and the chemical potential is given by Eq. (29) with the 
dimensionlcss parameter u obeying the equation 



(47) 



1/w + lnu = V2irl z /a — 7 — ln(167rn2D^ z )- 



This result for p is well consistent with relations (38) and 
(39). Indeed, substitution of Eq. (44) with constant (46) 
into Eq. (39) leads to Eq. (29) provided that the relation 
(37) is employed in conjunction with the approximation 
cxp[— {z\ — z 2 ) 2 /(2l 2 )] ~ 1 due to the integration with the 
short-range potential with R c <ti l z . In the leading order 
at small 2D densities, expressions (29) and (47) result in 



/! 



47rfi 2 7l2D 



1 



m \/2Trl z /a-~f -ln(16irn 2 T)l 2 ) 



(48) 



This differs from the result [35] of Petrov, Holzmann, 
and Shlyapnikov only by the additional numerical term 
—7— In 2 = —1.2703 ... in the denominator. We note that 
the healing length in two dimensions takes the form £ = 
l/^AimzDU, which differs from that in three dimensions 
l/\/4-7rna. Due to the criterion l/y^rTiB <C £, the 



obtained results relate to sufficiently small densities, for 
which uCl. 



3. ID regime 

Contrary to the 3D and 2D non-ideal Bose gases, there 
is no Bose-Einstein condensate in one dimension [10, 59] 
in the thermodynamic limit, because the long- wave fluc- 
tuations of the phase break the off-diagonal long-range 
order. Nevertheless, one can speak about the quasi- 
condensate [44] if a size of the of the ID system is suf- 
ficiently small. Indeed, at zero temperature the phase 
fluctuations are suppressed if ln(£ z /£) <C «id£ [44, 45], 
which can be fulfilled only at finite number of particles. 
Here L z stands for the size in z-direction. 

All calculations concerning the ID quasi-condensate 
in the case R e -C l p -C £ can be done in complete 
analogy with the 2D inhomogeneous Bose gas consid- 
ered in the previous subsection. The gas is strongly 
confined in the x-y plane by the harmonic trapping 



potential in V e > 



p j2 with the length l p 



y/H/(mojp), and remains homogeneous in z-dircction. In 
the regime involved, we can put <f>(p) = y / nm<^o (p) , 
(j>o(p) = cxp[— p 2 /(2l 2 )]/l p y/n is the ground state solution 
of the one-particle Schrodinger equation with the energy 
Eq = Hujp. Reasoning by analogy with Sec. IIIB2, we 
obtain 

M = »id / dzU(z), (49) 

where we introduce the even function U(z) = U(—z) 

U(z) = Jdp 1 dp 2 V(V(p 1 -p 2 ) 2 + z 2 ) 
x<p(pi,P2,z)<fo{pi)<f>o{p2)/niD- 
The ID analogue of Eq. (40) is the equation 



2m dz 2 



1>{z) = U(z) 



(50) 



for the function 

i>(z) = / dpidp 2 ip(pi, P2, z)<f> (pi)(f>o(P2)/ni-D. 

Equation (50) can be rewritten in the Lippmann- 
Schwinger form at /i — > (see discussion in Sec. Ill A) 

<p{z) = 1 + (m/h 2 ) fdz'U(z')\z-z'\/2 (51) 



for the function ip(z), defining as 

<f(z) = / dp 1 dp 2 ip(p 1 ,p 2 ,z)4> (p 1 )4> (p 2 )/n 1D . 



(52) 



Equations (49) and (51) give the asymptotics for _R C <C 
z<C 



(p(z) ~ 1 + mp\z\/(2niY)H 



(53) 



On the other hand, in the region r <C l p <C ^ we can use 
the analogue of Eq. (44) 



J") 



V(pi,P2,z) = Cip ( 3 ^(r)n 1I )(l)o(pi)(l)o(P2), 



(54) 



which leads to the asymptotics after the integration in 
Eq. (52) 



lp(z) c-C- C(y/7/2 - \z\/l)a/l p . 
Comparing Eqs. (53) and (55) yields 

C = 1/(1 - y/^/2a/l p ), 
2H 2 niY) a 1 



/i 



III 



(55) 

(56) 

(57) 



which differs from Olshanii's result [47] through the nu- 
merical factors \pH — 1.772... in the denominator in- 
stead of the constant 1.4603 . . . introduced by him. We 
note that in the paper [47] oj^ = yj2h/(mujp) — v2 l p in 
our notation. One can see that the criteria of applicabil- 
ity of the obtained results l p <§C £ and 1/niD <^ £ impose 
the following restriction on the ID density 



a 1 

y2 < "ID < - 



a 



(58) 



since ^ ~ l p j \f2anra in one dimension. 



IV. THE KINETIC, INTERACTION, AND 

EXTERNAL FIELD ENERGY OF THE 

TRAPPED BOSE GAS 

The simplest method for obtaining the values of the 
interaction energy (3), the kinetic energy (4), and the 
energy of interaction with an external field (5), is to apply 
the variational theorem. The latter can be formulated in 
general as follows. If a function fix) obeys the functional 
equation 



6F[{f(x)},\}/5f(x) = Q 



(59) 



with the functional F depending on the function f(x) 
and the parameter A, then the solution of Eq. (59) 
f(x) = fo(x, A) is also dependent on A. Nevertheless, 
when calculating the derivative of the stationary value of 
the functional with respect to A, we can take into consid- 
eration only the explicit dependence on this parameter 

dF[{f (x, A)}, A]/dA = dF[{f Q (x, A)}, \}/d\. (60) 
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This is obvious due to Eq. (59). 

The variational theorem (60) is still valid if the func- 
tional contains two or more functions. In our case, the 
functions can be associated with 4>{x\) and ip(x\, x<i) in- 
volved in the energy functional (15). Considering M as 
the parameter of the variational theorem, we come to 
the standard thermodynamic relation dE/dN = \j! = 
fi + Eq. One can rewrite this derivative in terms of 
the energy per particle e = E/N and the density of 
particles dE/dN = d(en)/dn, which gives the relation 
e = (1/n) f Q n dn' [i(n') + Eq. Then relations (29) and 
(47) lead to 



£2D — 2irh ri2T)u/m ■ 



h 2 



2ml 2 , 



(61) 



with u given by Eq. (47). In the same manner, we obtain 
from Eq. (57) [60] 



£lD 



H 2 niu a 




(62) 



Equations (61) and (62) give us the equilibrium value 
of the energy (15) per particle in the 2D and ID cases, 
respectively. In order to calculate the interaction energy 
with the help of the variational theorem, one can replace 
V — > XV and differentiate e with respect to A at A = 1. 
All we need to know is the derivative of the 3D scattering 
length, which reads [16, 61] 



A 



da 

'dX 



da 



dm 4ttH 2 



d 3 r[A)] 2 \V(r). (63) 



It is convenient to introduce one more characteristic 
length [16], the positive parameter b, 



h \ da 

b = a - x dx 



A=l 



So, we have 



£2Dint 



£lDint 



= -l/dV|v4 >)| 2 . 



i- b - 

a 



2ttH 7T-2D 9 V 2?r I 



H 2 niu a 



I 2 



a 

l- b - 
a 



(64) 
(65) 



u 2 in 



where we use the approximation u 2 /(l — u) 

Eq. (64) and restrict ourselves by the leading order in 

Eq. (65). 

With the same method, replacing V cx t — > XV c ^t (which 
is equivalent to I — ► 1/ vX) and differentiating, we arrive 
at the external energy per particle 



&2Doxt 



£lDoxt 



2irH 2 ri2D u 2 ( \/2ttI 



m 4 \ a 
H 2 niD a H 2 



-2 



h 2 



Ami 2 



2m I 2 



2ml 2 ' 



, (66) 



(67) 



In the same 


manner, we have 


^kin — 


—mc 


)e/c 


?m, 


which 


leads to 










£2Dkin — 


2irH 2 n2B 2irH 2 ri2i) 2 

u u 

m m 


EC 


Fhs 


h 


■') 


a 






+ V2^l z / b\ 
a \ a J 


h 2 






(68) 




Ami 2 ' 


£lDkin — 


h 2 n\D b h 2 ni]j 
m I 2 2m 


a 

l 2 + 2 


a 2 
ml 2 






(69) 



One can see that sum of the kinetic, external and inter- 
action energies equals to the total energy, as it should be. 
Note that the developed formalism allows us to calculate 
the interaction energy directly, starting from the expres- 
sion (11) and using Eq. (37), since we have the analytic 
expressions (44), (46), (54), and (56) for the short-range 
behaviour of the anomalous average. 

We note that the ratio b/a need not be small. In par- 
ticular, it is of order of ten for the realistic interaction 
potentials of alkali atoms [48]. We stress that the term 
with the length b appears in the mean interaction energy 
by virtue of the the short-range two-body correlations at 
the distances of order of a and in the mean kinetic en- 
ergy by sufficiently large momenta of order of p > H/a 
in the momentum distribution. In the static structure 
factor, this region is rather difficult to be measured ex- 
perimentally. The density expansion method gives the 
value of the release energy that is defined as sum of the 
interaction and kinetic energies 



£2Drol 



27rfi 2 7i2D 27rfi 2 rj 2 D u 2 ( s/2ttI 



£lDrol 



m 

h 2 
+ 4ml 2 z ' 
h 2 n\£) a 
2m 21 2 



2ml 2 ' 



(70) 



(71) 



As one can see, the parameter b is canceled and not in- 
volved in the values of the release energy. Let us compare 
the values of the release (70-71) and total energy (61-62). 
The energy of zero-point oscillation is involved in the re- 
lease energy with the factor 1/2, as it should be for the 
harmonic trap. The other terms would coincide in the 
standard GP approach, but we have obvious difference 
due to accounting for the non-condensate contribution. 
In principle, the obtained corrections should be measur- 
able in experiments. 



V. VIRIAL THEOREM 

The virial theorem can be obtained immediately from 
the energy functional (15) if we consider its variation 
in vicinity of the stationary state (ground state) with 
respect to the scaling transformation of the ground 
state functions 0o and ipQ, obeying the generalized GP 
equations (16) and (17). Namely, we substitute into 
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Eq. (15) the functions (j)( r i) = ct 3 ^ 2 4>o(ari) and VK r i) = 
a 3- 0o(o ;r ij otTi). Replacing the variables in the integrals 
ri — > ari and r2 — ► OT2, we notice that, first, the last 
term in the functional equals to zero for any a, and, 
second, the other terms can be written in terms of its 
stationary values 



E{a) = a 2 E kin + E^ t /a 2 + E int [V(r/a)] . 



(72) 



Since the variation of the functional should be zero for 
any small variations of the functions, we have dE/da = 
at a = 1, which leads to 

2£ kin - 2E cxt + E int [-rV'(r)} = 0, (73) 

where the terms are given by Eqs. (11)-(13). The value of 
the last term corresponds to the interaction energy with 
the potential —rV'(r) = — rdV(r)/dr. In the case of the 
GP approximation (36), one can simplify the last item in 
Eq. (73) by means of Eq. (37) and relation [6] 



47rfi 2 



in 



dr 4nr 2 fo>(°) (r)] 2 2V(r) + r^-^- 
V dr 



The result takes a form 



E int ~ iydR|^(R)| 4 |dr[-ry'(r)][^(r)] 2 

-(3a -26) /dR|0(R)| 4 . (74) 



2irH 2 



If the potential is of the weak-coupling type [18], one can 
neglect b <C a and arrive at the virial theorem obtained 
for the ^-function interaction potential [3] . 

If the system is homogeneous in the x-y plane (the 2D 
Bose gas of Sec. Ill B 2) or in the z direction (the ID Bose 
gas of Sec. IIIB3), it can be considered as confined by 
infinite walls in associated directions. Then one should 
be careful when deriving the virial theorem from Eq. (72), 
as all its terms relate to the density ri2T>/or or riiu/cx for 
the 2D or ID Bose gas, respectively. For this reason, we 
come to 

2n 2 D- = 2£ 2 Dkin- 2e 2 Doxt+ £2Dint [-rV (r)] , (75) 



"ID 



dn 2 D 
de 1D 

dn 



2e 



lDkii 



ID 



2eiDext+eiDint[-rV' / (r)].(76) 



The interaction term in these equations can be easily cal- 
culated by analogy with Eq. (74) but using Eqs. (44) and 
(54), respectively. It is not difficult to be convinced with 
the help of Eqs. (61) and (62) that the virial theorems 
(75) and (76) are fulfilled. 

One can also find a relation between the chemical po- 
tential and the various parts of the energy. Let us multi- 
ply Eq. (16) by 4>{x\) and integrate over x\, and multiply 
Eq. (17) by ip(x±,X2) and also integrate over x\ and Xi- 
Summing the obtained expressions yields 

N/J, = -Ekinl + 2i?kin2 + £"extl + 2i? oxt 2 + 2_E; nt , (77) 



Here, E oxt i and -Ekini are the condensate contributions 
in the external and kinetic energies given by the last 
terms in Eqs. (12) and (13), respectively, and E cxt 2 and 
Ekin2 are associated with the non-condensate contribu- 
tions, given by the residual parts of these equations. One 
can easily see that the relation (77) is fulfilled with _E C xti 
and -Ekini corresponding to the last terms in Eqs. (66) 
and (67), and (68) and (69) for the 2D and ID Bose 
gases, respectively. One can notice that £^112 could be 
negative for the ID Bose gas, if b < a/2 [see the first two 
terms in Eq. (69)]. Certainly, this is not a drawback of 
Eqs. (16) and (17) it is but due to the choice of anzatz 
4>{p) = \ /n i~D ( t'o(p)i which leads to overestimation of the 
quasicondcnsate contribution -Ekini in the ID kinetic en- 
ergy. Indeed, the Gaussian profile ™id|</>o(p)| 2 relates to 
the total density of the ID gas ($t(r)&(r)) but not to 
the "quasicondensate component" ^(p)! 2 . The latter is 
difficult to define accurately in the ID case, since there is 
no eigenvalue of the one-body density matrix that is pro- 
portional to the total number of particles. Nevertheless, 
we stress that the total value of -EiDkin is positive, and 
the results (65), (67), and (69) look quite reasonable. 



VI. CONCLUSIONS 

The main result of this paper are the generalized 
GP equations in the time-dependent (18-19) and sta- 
tionary form (16-17), which allow us to determine the 
interaction term self-consistently for interaction poten- 
tials even containing a hard-core. The method, which 
can be used for homogeneous, strongly inhomogeneous 
quasi-low-dimcnsional, and cross-over regimes was de- 
rived within a general HFB framework. 

The HFB method is a mean-field approximation, 
which generally works well only for weak-coupling po- 
tentials [18]. In order to extend the HFB scheme 
to hard-core potentials, the bare interaction potential 
is usually replaced by a renormalized pseudopotential 
V(r) —* (47rfi 2 /m)(5 3 (r). However, such a replacement 
leads to an ultraviolet divergence and incorrect treat- 
ment of short-range correlations of the particles. We 
have shown that the appropriate renormalization can be 
obtained within the HFB scheme if, from the two-body 
density matrix, only the anomalous correlation function 
ip(xi,X2) = (v& (xi)$ (X2)) is retained. The anomalous 
correlation function can be interpreted as the wavefunc- 
tion of two bosons in the condensate. Its short-range 
behaviour is described well in the proposed scheme at 
the cost of loosing the correct description of the long- 
range behaviour. However, long-range correlations are 
not needed for deriving the non-linear term in the gen- 
eralized GP approach, which instead is determined by 
short-range correlations. Methods which can describe 
both the short- and long-range correlations accurately 
were discussed in Refs. [13, 16, 17, 24], but these meth- 
ods are appropriate only for the homogeneous Bose gas. 
The method proposed in this paper was shown to work as 
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well in inhomogeneous situations. Cigar (quasi-lD) and 
pancake (quasi- 2D) geometries were considered as exam- 
ples. Furthermore, it was shown that the contribution 
of short-range correlations to the kinetic and release en- 
ergies of a tightly trapped gas can be calculated within 
this scheme and that they are substantial. Interesting 
future applications of the proposed method may include 
the modification of the nonlinearity in quasi- ID waveg- 
uides [62, 63] and molecular Bose condensates in optical 
lattices. 
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APPENDIX A: TWO-BODY WAVE FUNCTIONS 

IN THE HARTREE-FOCK-BOGOLIUBOV 

APPROXIMATION 

In general, the two-body density matrix can be ex- 
panded in a complete set of its eigenfunctions 



(^(x 1 )^(x 2 )^(x , 2 mx , 1 )) = J2 N ^<M-* 



x 2 ) 



XLp u ^(x[,x 2 ). (Al) 

The eigenfunctions can be called two-body or pair wave 
functions. If they are normalized to unity, it follows from 
Eq. (Al) that fdx 1 dx 2 (&(x 1 )&(x 2 )i'(x 2 )i>(x 1 )) = 
N(N - 1) = ^2 N v>li , i.e., the sum of all N v ^ is the 
total number of pairs. Therefore, the non-negative quan- 
tity N v ^ can be interpreted as the mean number of the 
pairs in the state (v, /i) , any pair being doubly taken. 

Let us consider the homogeneous spinless Bose gas in 
the HFB approximation [22, 25]. Within that approxi- 
mation, the two-body wave functions can be easily cal- 
culated [15]. The statistical average of any product of 
quantum operators i? and ffi can be calculated with the 
Wick-Bloch-Dc Dominicis theorem [64] , since the Hamil- 
tonian is approximated by a quadratic form of the Bose 
operators & P and a p connected with initial operators a p 
and a p by the canonical Bogoliubov transformations (see 
Appendix B). Extracting the c-number part \& = ^/riQ+'d 
and ffi = ^/tlq + ffi and using that theorem, one can 
rewrite the four-boson average in the form 



*.(r + £)*.(r-E)*(r<-1)*(r' + 1) 

2- «(q/2 - p)^j 

v2cos(p • r)v2 cos(p • r') 
(A2) 



= nl<p*(r)if(r') + / d 3 pd 3 q 

n(q/2 + p) n(q/2 - p 
(2tt) 3 (2tt) 3 

x exp[zq • (R' — R)], 



where we put by definition <f>(r) — 1 + (^(R + r/2)$(R — 
r/2))/no- Because the expansion (A2) is written in the 
thermodynamic limit, the sum in Eq. (Al) becomes an 
integral. The Bose-Einstein condensate manifests itself 
in presence of <5-functions in this integral (note that the 
first term in the r.h.s. can be included in the integral 
with the help of the 5- functions). By comparing the rep- 
resentation (Al) with that of (A2), one can conclude the 
following: 

(i) The quantum numbers of the pair wave functions are 
the relative momentum v = p and the center-of-mass 
(total) momentum /i = q of two particles; all these func- 
tions belong to continuous spectrum and thus describe 
the scattering of two bosons in the medium of the other 
bosons. 

(ii) The maximum eigenvalue Nq(Nq — 1) — Nq with 
p = q = corresponds to the state of two par- 
ticles in the condensate; its normalized eigenfunction 
ip(r)/V can be interpreted as a pair wave function of 
the condensate-condensate type. Thus, the anomalous 
average (#(r)$(0)) can be associated with the scatter- 
ing part of the two-body wave function of the bosons in 
the condensate [15]; in particular, it is responsible for 
the short-range spatial correlations of two bosons in the 
Bose-Einstein condensate. 

(Hi) The other macroscopic eigenvalues 2Non q with q = 
±2p correspond to the two-body states with one particle 
in the condensate and another one beyond the conden- 
sate; its eigenfunctions v / 2cos(q-r/2) exp[iq-R]/V are of 
the condensate-noncondensate type [66]. The residuary 
non-macroscopic eigenvalues n(q/2 + p)n(q/2 — p) are 
related to the noncondensate-noncondensate pairs with 
the two-body wave functions y/2 cos(p • r) exp[iq • R]/V. 



Note that the wave function of the condensate- 
condensate type is not reduced to a product of two one- 
body wave functions in the condensate, which equal to 
1/vV for the homogeneous Bose gas. This is obvious, as 
particles in the Bose-Einstein condensate interact with 
each other and with the other particles beyond the con- 
densate. Another important point is that all the other 
two-body wave functions are symmetrized plane waves 
(consistent with the Born approximation) in the frame- 
work of the HFB method. This is evidently a disadvan- 
tage of the HFB scheme. As a consequence, we always 
arrive at divergences for a hard-core potential when eval- 
uating the contribution of the condensate-noncondensate 
and noncondensate-noncondensate wave functions in the 
interaction energy (3). At the same time, the contribu- 
tion of the condensate-condensate "channel" should be 
finite in the interaction energy provided the anomalous 
averages are calculated in a self-consistent manner. The 
generalization of the expansion (A2) beyond the HFB 
approach and more detailed discussions can be found in 
Ref. [15]. The pair wave function method of Ref. [15] was 
generalized to the inhomogeneous systems in Ref. [68]. 
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APPENDIX B: RELATION BETWEEN THE 

NORMAL AND ANOMALOUS TWO-BOSON 

AVERAGES 

Let us establish a relation between the normal 
(w (xi)i!}(x2)} and the anomalous average (i}(xi)'&(x2)) 
for the vacuum state, which describes the behaviour of 
the N-body system at zero temperature, in the frame- 
work of the Hartree-Fock-Bogoliubov method. We re- 
member that the vacuum state |0) is defined as a„|0) =0 
for any v ^ 0, here the quasiparticle creation and de- 
struction operators a\, and d„ can be introduced through 
the Bogoliubov transformation (/ ^ 0) 



d/ = J^ {uf v a v + Vf v a\,), 

V 

4 = E ( u }» ,i t+ v % & v)> 



(Bl) 

(B2) 



where / and v denote discrete (multi-)indices. The sum 
^2 v means ^2 u u. and the Bose-operators d , and d/ cre- 
ate and destruct a particle in the eigenstate </>f (x) of the 
one-body matrix (^>'(x')^(x)) 



dx' (W(x')${x))<l>f(x') = n f <t> s (x), 



normalized as J dx \4>f(x)\ 2 = 1. Note that the set of 
cigenfunctions including the normalized condensate func- 
tion 4>o(x) = ($ (x)} J \fNo with A^o = n/=o is complete 
and orthogonal 



E ( t ) *f( x ) < t ) f( x ') = 5 ( x ~ x ')> 
where we define the "discrete" ((-function as 

w-{S:/7!!: 



(B3) 

(B4) 



S(f 



From the Bosc commutation relations [d/ , df , 

f) and [af,a\f,\ = S(f-f') and Eqs. (B1-B2) we obtain 

2'(«/ v «>v-v/^K) = *(/"/'). (B5) 
2^ {u fv Vf v -Vf v Uf> v ) = 0. (B6) 

By using the definition of the quasiparticle vacuum state 
and Eqs. (B1-B2), we can calculate the averages 

F(f,f) = {a)a f ,)=Y^ v} v v f . v , (B7) 

V 

$(/,/') = (a/a,,) =£)' «,„«/,„. ^ B8 ) 



Our purpose is to find the relation between the normal 
F(f, /') and the anomalous $(/, /') averages for that 
state. In order to simplify our calculations, we rewrite 
Eqs. (B1-B2) in the matrix notations 



X 



, X 



U V 

v* u* 



(B9) 



Here the matrix X is composed of the matrix (U)i 



and (V)i 



The columns contain the operators 



d/ and dl, and a u and aj,, respectively. We use the 



standard notations for the complex conjugate {V*)ij 
v* ■ , transposed (V 1 -)^ — Vji, and Hcrmitian conjug. 
matrix (V*)^ = v* % . Then Eqs. (B7-B8) read 

F = V*V T = F\ $ = UV t = <P 
and Eqs. (B5-B6) can be written as 

u v\(W ~v T 
v* u* l-vt U T 



(BIO) 



1 
1 



(BH) 



where 1 denotes the identity matrix. Let us introduce 
the composed matrices 



CT3 
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-1) ' a+ 10 



1 



and rewrite Eq. (Bll) in the form 
Xa 3 X^a 3 = 1, 



(B12) 



(B13) 



where 1 stands now for the composed identity matrix, i.e. 
the r.h.s. of Eq. (Bll). The matrix representation (B13) 
is very convenient. For example, from this equation we 
have X^ 1 = CT3XV3, and 



a 3 XV 3 L a t =(_ Ft uT 



which reads in usual notations 



a f 



V 






E 



{u* vf a v - Vufal), 
{u uf al - v* vf a v ). 



This equation together with the commutation relations 
leads to 

2^ (KfUuf - v v fv* vfl ) = 6(f-f), 



E 



which is nothing else but the matrix equation 
X^a 3 Xa 3 = 1, resulting from Eq. (B13). 

Employing the idea of Ref . [65] , in which the Hartree- 
Fock-Bogoliubov method for Fermi systems was devel- 
oped, we define the matrix K with the help of the nota- 
tions (BIO) and (B12) 



K = X f a 3 a + X<T 3 



t + F* -$ 



([)* 



-F 
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Due to Eq. (B13) and the relation (cr + ) 2 = c + we have 
K 2 = K . Rewriting the latter equation in terms of the 
matrix F and $, we obtain two independent relations 



$*$ = F + F z 



(B14) 
(B15) 



F*$ = $F, 
which read in components 

E^* (A ./)*(/- h) = ^V(/ 1 ,/)F(/,/ 2 )+F(/ 1 ,/ 2 ), 

(B16) 

^V(/, /!)$(/, / 2 ) = £$(fi,f)F{f,f 2 ). (B17) 
/ / 

By using these equations, Eqs. (B3-B4), and the defini- 
tion $(x) = Y] ,, a v <f> v {x), one can rewrite Eqs. (B16- 
B17) in the coordinate representation 

dx {$(xi)#(x))(ft(x)'&(x2)), (B18) 



dx(^(x)d(x 1 ))(^(x)d(x 2 )) = / dx (0(xi)0(x)) 

x (&(x)S(x 2 )). (B19) 



If the condensate depletion is small, one can neglect the 
second term in the r.h.s of Eq. (B18), which is of the next 
order. Thus, we obtain the expression 



(&( Xl )S(x 2 )) ~ dx(^(x 1 )^(x)}0(x)^(x 2 )}. (B20) 



Note that Eq. (B19) turns into identity in the approxi- 
mation (B20), and the same is valid for Eqs. (B15) and 
(B17). 
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